A Multi-Model Based Stability Analysis Employing Multi-Environmental Trials (METs) Data for Discerning Heat Tolerance in Chickpea (Cicer arietinum L.) Landraces

Identifying a congenially targeted production environment and understanding the effects of genotype by environmental interactions on the adaption of chickpea genotypes is essential for achieving an optimal yield stability. Different models like additive main effect and multiplicative interactions (AMMI 1, AMM2), weighted average absolute scores of BLUPs (WAASB), and genotype plus genotype–environment (GGE) interactions were used to understand their suitability in the precise estimation of variance and their interaction. Our experiment used genotypes that represent the West Asia–North Africa (WANA) region. This trial involved two different sowing dates, two distinct seasons, and three different locations, resulting in a total of 12 environments. Genotype IG 5871(G1) showed a lower heat susceptibility index (HSI) across environments under study. The first four interactions principal component axis (IPCA) explain 93.2% of variations with significant genotype–environment interactions. Considering the AMMI stability value (ASV), the genotypes IG5862(G7), IG5861(G6), ILC239(G40), IG6002(G26), and ILC1932(G39), showing ASV scores of 1.66, 1.80, 2.20, 2.60, and 2.84, respectively, were ranked as the most stable and are comparable to the weighted average absolute scores of BLUPs (WAASB) ranking of genotypes. The which–won–where pattern of genotype plus genotype–environment (GGE) interactions suggested that the target environment consists of one mega environment. IG5866(G10), IG5865(G9), IG5884(G14), and IG5862(G7) displayed higher stability, as they were nearer to the origin. The genotypes that exhibited a superior performance in the tested environments can serve as ideal parental lines for heat-stress tolerance breeding programs. The weighted average absolute scores of BLUPs (WAASB) serve as an ideal tool to discern the variations and identify the stable genotype among all methods.


Introduction
Chickpea (Cicer arietinum L.) is an important grain legume crop with a high protein content, implying its significance in sustaining food and nutritional security for inhabitants of arid and semiarid regions [1][2][3].Chickpea-growing regions are predominantly located in arid and semiarid areas, making them susceptible to high-temperature stress, mainly during the reproductive and grain-filling stages [4][5][6][7][8].The impact of heat stress on chickpea yield and production has become a matter of significant concern in recent times.The paradigm shift in chickpea cultivation is driving factors including a shift in major cropping areas from cooler long-season climates to shorter warm-season climates, the expansion of late-sown chickpea areas due to higher cropping intensities, and also the overall rise in seasonal temperatures caused by global climate change [9][10][11].From germination to the vegetative and reproductive phases, all stages of chickpea growth are hindered by temperatures exceeding the normal range [12][13][14][15][16][17].During the terminal stage of chickpea growth, which encompasses the reproductive phase, heat stress can have significant consequences for both crop yield and quality.This critical stage involves the formation of flowers and pods and the development of seeds [18][19][20].
There is an urgent need to identify stable sources of donors that can be utilized into a breeding program [21,22].The identification of landraces known for their acclimatization under harsh and hardy environments is a basic prerequisite.Testing their response to Indian agroclimatic conditions for normal and late-sown conditions can provide crucial clues about genotypes with inherent heat-stress-tolerant capacities [23,24].Genotypeenvironment interaction (GEI) occurs when different genotypes show varying responses to changing environments with fluctuations in crop yield.Researchers rely on stable statistical measures that include both univariate and multivariate techniques.These methods encompass cluster analysis, pattern analysis, and principal component analysis (PCA) or biplots.Biplots, in particular, are highly effective graphical tools used to represent the relationships between genotypes and environments visually to analyze genotype-environment interactions (GEIs).Two widely utilized biplot models are the AMMI biplot (additive main effects and multiplicative interaction) and the GGE biplot (genotype-plus-genotype × environment) [25][26][27][28][29]. Researchers commonly employ genotype-plus-genotype by environment (GGE) biplots for various purposes [30][31][32], such as classifying mega environments, assessing genotype rankings, and selecting discriminative and representative environments.Best linear unbiased prediction (BLUP) had a higher predictive accuracy than any AMMI family member in analyzing multi-environment trials (METs) with a random analysis of the genotype-environment interaction (GEI) effect using a linear mixed-effect model (LMM).The weighted average absolute scores of BLUPs (WAASB) can reliably estimate stability in a bidimensional plot with multiple interactions on the principal component axes (IPCAs), which is essential for ranking genotypes.The simultaneous selection index WAASBY is useful when different weights are needed for performance and stability [33][34][35][36].
In this study, we used multi model analysis to identify high-yielding, stable chickpea genotypes across diverse environments using various stability analyses.

Results
According to Fischer and Maurer (1978) [37], the HSI (heat susceptibility index) has served as the predominant criterion for identifying tolerant genotypes.In our study, HSI was calculated using timely and late-sown data of three locations for two consecutive years.Results showed that IG 5871 (0.12, 0.12, 0.18, 0.07, 0.1, 0.08) had a lower HSI across environments and across locations, along with ILC 8666, IG 5905, and IG 5868, which exhibited lower yield dips under more than four environments (Supplementary Table S1).

Variance Analysis for Yields
The combined analysis of variance (Table 1) revealed highly significant differences (p ≤ 0.001, p ≤ 0.01) for yield per plot among different factors: locations, genotypes, and the genotype-location (G × E) interaction.Yield variations result from evolving environments and genetic diversity.Each violin plot (Figure 1) consists of two mirrored distributions, one for normal-sown conditions (blue color) and the other (green color) for late-sown conditions.Yield variations result from evolving environments and genetic diversity.Each violin plot (Figure 1) consists of two mirrored distributions, one for normal-sown conditions (blue color) and the other (green color) for late-sown conditions.The width of each violin corresponds to the density of data points, with broader sections indicating higher data density and narrower sections indicating lower density.The plots illuminate differences in central tendencies, spread, and distribution shapes, offering valuable insights into the impact of late conditions on crop yields.Genotype performance variation was prominent among locations, with 18.47% attributed to location, 4.03% to replication, 16.19% to genotype, and the highest, 21.30%, to genotype-environment interaction across locations.Moreover, qualitative (crossover type) interaction was demonstrated where the rank order of genotypes changed across environments (Supplementary Figure S1A and Table S2), indicating the presence of substantial genotype-environment The width of each violin corresponds to the density of data points, with broader sections indicating higher data density and narrower sections indicating lower density.The plots illuminate differences in central tendencies, spread, and distribution shapes, offering valuable insights into the impact of late conditions on crop yields.Genotype performance variation was prominent among locations, with 18.47% attributed to location, 4.03% to replication, 16.19% to genotype, and the highest, 21.30%, to genotype-environment interaction across locations.Moreover, qualitative (crossover type) interaction was demonstrated where the rank order of genotypes changed across environments (Supplementary Figure S1A and Table S2), indicating the presence of substantial genotype-environment interactions.The overall average yield recorded 174.5 g/plot, with the minimum value of 110 g/plot in the 2021 Dharwad late-sown condition, and the maximum yield value of 257 g/plot recorded in the 2022 Delhi timely-sown condition (Supplementary Figure S1B).Annichiarico environmental index was used to predict favorable and unfavorable environments along with the most suitable genotype for the particular environment (Supplementary Table S3).Descriptive statistics, variance components, and genetic parameters for grain yield of all the genotypes evaluated across the year and locations are provided in Supplementary Table S4, where overall heritability (broad sense) across environments was reported to be 72.54%.To analyze genotype-environment interaction (GEI), a principal component analysis (PCA) was executed (Table 1).The ordination method using an approximate F-statistic showcased significant distinctions across PC1 (35.4%),PC2 (28.9%), and PC3 (21.7%), totaling 86% variation.Inclusion of PC4 (7.2%) led to a cumulative 93.2% explanation of variability, with a significant (p < 0.001) effect.

Additive Main Effects and Multiplicative Interaction: AMMI 1
The biplot displays the first principal component (PC1) term on the abscissa and significant trait influence on the ordinate.Environments E3, E9, and E10, representing plot yield, showed PCA1 scores closer to zero, parallel to the average yield line.This suggests the strong performance of all genotypes in these environments.E6, E5, and E1 had vectors parallel to PC1 (35.4%), indicating a higher contribution to overall variation.E4, E2, E8, E1, E7, E9, and E3 had higher average yields than E10, aligned with the average yield line. E12, E6, E5, and E11 were below-average environments (Figure 2A).Certain genotypes like IG5862(G7) and ILC1312(G36) consistently excelled across all environments, with high mean yields near the axis origin.Genotypes ILC1312(G36), IG5856(G4), IG5904(G18), and IG5980(G21) also had higher mean yields, while IG5858(G5) and ILC1313(G37) had lower yields, aligning with zero scores on the first PCA1 axis.These genotypes were less affected by environmental variations.
interactions.The overall average yield recorded 174.5 g/plot, with the minimum value of 110 g/plot in the 2021 Dharwad late-sown condition, and the maximum yield value of 257 g/plot recorded in the 2022 Delhi timely-sown condition (Supplementary Figure S1B).Annichiarico environmental index was used to predict favorable and unfavorable environments along with the most suitable genotype for the particular environment (Supplementary Table S3).Descriptive statistics, variance components, and genetic parameters for grain yield of all the genotypes evaluated across the year and locations are provided in Supplementary Table S4, where overall heritability (broad sense) across environments was reported to be 72.54%.To analyze genotype-environment interaction (GEI), a principal component analysis (PCA) was executed (Table 1).The ordination method using an approximate F-statistic showcased significant distinctions across PC1 (35.4%),PC2 (28.9%), and PC3 (21.7%), totaling 86% variation.Inclusion of PC4 (7.2%) led to a cumulative 93.2% explanation of variability, with a significant (p < 0.001) effect.

Additive Main Effects and Multiplicative Interaction: AMMI 1
The biplot displays the first principal component (PC1) term on the abscissa and significant trait influence on the ordinate.Environments E3, E9, and E10, representing plot yield, showed PCA1 scores closer to zero, parallel to the average yield line.This suggests the strong performance of all genotypes in these environments.E6, E5, and E1 had vectors parallel to PC1 (35.4%), indicating a higher contribution to overall variation.E4, E2, E8, E1, E7, E9, and E3 had higher average yields than E10, aligned with the average yield line. E12, E6, E5, and E11 were below-average environments (Figure 2A).Certain genotypes like IG5862(G7) and ILC1312(G36) consistently excelled across all environments, with high mean yields near the axis origin.Genotypes ILC1312(G36), IG5856(G4), IG5904(G18), and IG5980(G21) also had higher mean yields, while IG5858(G5) and ILC1313(G37) had lower yields, aligning with zero scores on the first PCA1 axis.These genotypes were less affected by environmental variations.

Figure 2. (A)
The "AMMI 1" biplot depicts the primary effects of traits, and the first principal component (PC1) impacts of both genotype and environment.(B) The "AMMI 2" biplot visualizes the combined effects of the first two principal components (PC1 and PC2) for genotypes along with the genotype-environment interaction impact of 42 genotypes across two seasons and three locations, considering yield per plot (g/plot).

Additive Main Effects and Multiplicative Interaction: AMMI 2 Figure 2. (A)
The "AMMI 1" biplot depicts the primary effects of traits, and the first principal component (PC1) impacts of both genotype and environment.(B) The "AMMI 2" biplot visualizes the combined effects of the first two principal components (PC1 and PC2) for genotypes along with the genotype-environment interaction impact of 42 genotypes across two seasons and three locations, considering yield per plot (g/plot).

Comparison of BLUP and AMMI Family Models
The evaluation determines the superior model in a given context.Based on our diverse datasets exhibiting various genotype-environment interaction (GEI) patterns, our analysis suggested that best linear unbiased prediction (BLUP) emerged as the most predictively precise model (Figure 3B).Furthermore, we observed that AMMI0 demonstrated the highest accuracy within the AMMI models.Among the best linear unbiased prediction (BLUP) models of environment, genotypes, and genotypes by environment, in conjunction with AMMI (0-10 and AMMIF), our observations point to the best linear unbiased prediction (BLUP) by genotypes as the most accurate model.

Comparison of BLUP and AMMI Family Models
The evaluation determines the superior model in a given context.Based on our diverse datasets exhibiting various genotype-environment interaction (GEI) patterns, our analysis suggested that best linear unbiased prediction (BLUP) emerged as the most predictively precise model (Figure 3B).Furthermore, we observed that AMMI0 demonstrated the highest accuracy within the AMMI models.Among the best linear unbiased prediction (BLUP) models of environment, genotypes, and genotypes by environment, in conjunction with AMMI (0-10 and AMMIF), our observations point to the best linear unbiased prediction (BLUP) by genotypes as the most accurate model.
most elevated predicted means (207.9, 233.7, 228.7, 200.1) (Figure 4), consequently earning the distinction of "universal winners" due to their minimal IPCA1 scores (−5.59, −4.77, −4.36, and −3.9) (Supplementary Table S6).When examining genotype ranking based on the number of retained interaction principal component axes (Figure 4A-C), eleven interaction principal component axes (IPCAs) were considered for the chickpea dataset (Supplementary Table S6).It is apparent that the genotype ranking was significantly impacted by the number of interaction principal component axes (IPCAs) employed in weighted average absolute scores of BLUPs (WAASB) estimation, particularly up to four interaction principal component axes (IPCAs).Clear groups of genotypes with akin stability performance are discernible by the color variations on the left side.For instance, genotypes ILC239(G40), IG6002(G26), IG5862(G7), ILC0(Czech)(G30), and IG5861(G6) showcase the lowest WAASB values (0.77, 0.93, 1.03, 1.10, and 1.21) when considering four or more interaction principal component axes (IPCAs) (Supplementary Table S6).They are consecutively ranked as the first, second, third, fourth, and fifth most stable.Parallel rankings are observed with the WAAS index for these genotypes (0.98, 1.27, 1.35, 1.46, and 1.54) in the first cluster (dark blue gradient) as we observed a correlation between ASV and weighted average absolute scores of BLUPs (WAASB) (Supplementary Figure S4).However, when accounting for ASV, genotypes IG5862(G7), IG5861(G6), ILC239(G40), IG6002(G26), and ILC1932(G39) receive ASV scores of 1.66, 1.80, 2.20, 2.60, and 2.84, respectively.We observed significant correlation among stability models (Supplementary Figure S4).Consequently, ILC1932(G39) would be ranked as the fifth most stable (Supplementary Figure S3), although in actuality, it belongs to the second cluster of genotypes (red color) with the highest stability.The "which-won-where" pattern of multi-environment trials (METs) vividly illustrates genotype-environment interaction (GEI) by correlating genotypes and environments.It also aids in evaluating genotypes stability versus mean performance across environments and assessing test environments for representativeness and discrimination  The "which-won-where" pattern of multi-environment trials (METs) vividly illustrates genotype-environment interaction (GEI) by correlating genotypes and environments.It also aids in evaluating genotypes stability versus mean performance across environments and assessing test environments for representativeness and discrimination (Figure 5A).The polygon view depicting the "which-won-where" pattern showcases the combined impact of genotype main effects, with the genotype ILC8666(G41) and IG5896(G17) positioned along the line connecting JG14(G42) and ILC1932(G39) affirmed the ranking order JG14(G42) > ILC8666(G41) > IG5896(G17) > ILC1932(G39) across all environments.These equality lines partitioned the biplot into sectors, with the leading genotype for each sector located at the corresponding vertex.In this instance, the twelve environments were grouped into one sector.JG14(G42) emerged as the universal winner, particularly in E9 and E10, while ILC1932(G39) exhibited subpar performance in these environments displaying the crossover (Figure 5A,B).The GGE biplot polygon view depicting the "which-won-where" pattern showcases the combined impact of genotype main effects and G × E interaction for 42 C. arietinum genotypes across two seasons and three locations in terms of yield per plot.These biplots were generated using centering = 2, SVP = 3, and scaling = 0 settings.

The Assessment of Stability across Environments Using the "Mean vs. Stability" Analysis of the GGE Biplot
The assessment of stability across environments using the "mean vs. stability" analysis of the genotype plus genotype by environment (GGE) biplot is an effective approach.In Figure 6A, the abscissa of the average-environment coordination (AEC) indicates a higher mean yield across various environments.Genotypes such as JG14 (G42), IG5871(G1), and IG5868(G11) exhibit the highest mean yields, while closely trailing are IG5895(G16), IG5865(G9), IG5852(G3), and others.IG5980(G21) demonstrates a mean yield similar to the grand mean, in contrast to ILC1932(G39), which exhibits the lowest mean yield.The AEC ordinate, represented by the vertical line, signifies greater variability (indicative of poorer stability) in either direction.Notably, IG5980(G21) and ILC1932(G39) are characterized by high instability, while IG5866(G10), IG5865(G9), IG5884(G14), and IG5862(G7) showcase considerable stability.The instability of ILC1932(G39) stems from its unexpectedly lower yield in several environments (E1, E7, E8, E2, E4, E9, E3, and E10) juxtaposed with higher-than-anticipated yield in E12, E11, E6, and E5.The GGE biplot polygon view depicting the "which-won-where" pattern showcases the combined impact of genotype main effects and G × E interaction for 42 C. arietinum genotypes across two seasons and three locations in terms of yield per plot.These biplots were generated using centering = 2, SVP = 3, and scaling = 0 settings.

The GGE Biplot s "Discriminativeness vs. Representativeness"
The genotype plus genotype by environment (GGE) biplot s "discriminativeness vs. representativeness" pattern serves as a crucial tool for discerning the most suitable test environments to select superior genotypes.The pattern discriminates among genotypes (discriminativeness), and, secondly, faithfully represents the overall set of evaluated
2.6.4.The GGE Biplot's "Discriminativeness vs. Representativeness" The genotype plus genotype by environment (GGE) biplot's "discriminativeness vs. representativeness" pattern serves as a crucial tool for discerning the most suitable test environments to select superior genotypes.The pattern discriminates among genotypes (discriminativeness), and, secondly, faithfully represents the overall set of evaluated environments (representativeness).In Figure 7A, our focus was on yield evaluations across distinct research locations (E9, E8, E2, E11, and E7).The concise vectors associated with these environments offer a clear insight into their discriminative potential.Environments characterized by expansive vectors, such as E5 and E6, exert a stronger discriminatory influence due to their larger variability.For instance, as depicted in the case of environment two (E10, E3, E12, E9) for yield, a subtle angle alongside an extended vector indicates an environment's heightened representativeness and discriminativeness.Moving to Figure 7A, environment rankings distinctly identify E10 and E3 as prime candidates for yield-focused genotype selection.Conversely, E5, E6, E7, and E1 are relegated to weaker positions, suggesting limited suitability for this purpose.In our evaluation of the three assessed locations, Dharwad emerges as the optimal choice, followed by Amlah and then Delhi.Our analysis of the angles between the AEC abscissa and the test environment vector further accentuates the differentiation between favorable and less ideal environments.Smaller angles, exemplified by environments like E10, E3, E9, E4, E12, E11, and E2, signify more advantageous conditions for genotype selection.Conversely, wider angles, observed in environments such as E5, E6, E1, E8, and E7, imply reduced potential for accurate genotype discrimination.Finally, the length of each environment vector provides a visual representation of its discriminatory capacity.Longer vectors, such as those associated with E5, E6, and E1, are more effective at distinguishing genotypes, emphasizing their significance in the selection process.Conversely, shorter vectors like E3, E10, and E12 display a diminished discriminatory ability.
Plants 2023, 12, x FOR PEER REVIEW 11 of 17 environments (representativeness).In Figure 7A, our focus was on yield evaluations across distinct research locations (E9, E8, E2, E11, and E7).The concise vectors associated with these environments offer a clear insight into their discriminative potential.Environments characterized by expansive vectors, such as E5 and E6, exert a stronger discriminatory influence due to their larger variability.For instance, as depicted in the case of environment two (E10, E3, E12, E9) for yield, a subtle angle alongside an extended vector indicates an environment s heightened representativeness and discriminativeness.Moving to Figure 7A, environment rankings distinctly identify E10 and E3 as prime candidates for yield-focused genotype selection.Conversely, E5, E6, E7, and E1 are relegated to weaker positions, suggesting limited suitability for this purpose.In our evaluation of the three assessed locations, Dharwad emerges as the optimal choice, followed by Amlah and then Delhi.Our analysis of the angles between the AEC abscissa and the test environment vector further accentuates the differentiation between favorable and less ideal environments.Smaller angles, exemplified by environments like E10, E3, E9, E4, E12, E11, and E2, signify more advantageous conditions for genotype selection.Conversely, wider angles, observed in environments such as E5, E6, E1, E8, and E7, imply reduced potential for accurate genotype discrimination.Finally, the length of each environment vector provides a visual representation of its discriminatory capacity.Longer vectors, such as those associated with E5, E6, and E1, are more effective at distinguishing genotypes, emphasizing their significance in the selection process.Conversely, shorter vectors like E3, E10, and E12 display a diminished discriminatory ability.

Test Environment Discrimination and Representation
Figure 7B illustrates distances among test environments for all traits.Concentric circles in the biplot visualize environment vector lengths, proportional to standard deviations within environments, reflecting environment discrimination capability.E5 and E6 are most informative, whereas E10 and E9 are least discriminative.Representativeness of test environments: E10 is highly representative, unlike E5 and E6 which are less so.Environments that are both discriminative and representative, like E10, suit broad genotype selection.Nonrepresentative yet discriminative test environments (e.g., E5, E6) aid specific megaenvironment genotype selection.From Figure 7B rankings, preferred target environments are E10 > E3 > E4 > E12 > E9 > E11 > E2 > E8 > E1 > E7 > E6 > E5.When focusing on a single mega environment, discriminative but nonrepresentative test environments (e.g., E5) effectively eliminate unstable genotypes.With regard to optimal environments for selecting broadly adapted genotypes, in a singular mega environment context, E10 emerges as the premier choice due to its proximity, while E5 and E6 prove suboptimal for region-wide cultivar selection.Notably, E5 and E6, sourced from Delhi_2021 timelyand late-sown seasons, contrast with E10 from Dharwad_2022 late-sown environment.Interrelationships among test environments (Supplementary Figure S5B) portray connections between test environments.E10, E3, E9, E2, E8, and E1 exhibit positive correlation (acute angle), unlike E5 and E1, with no apparent correlation (right angle).The distance between environments quantifies their distinctiveness in genotype differentiation.Nine environments cluster into two groups: E12, E11, E6, and E5 in one, and the remaining in another.Covariance, influenced by vector length and angle cosine, defines similarity between environments.Genotype assessment and ranking (Supplementary Figure S6A) utilizes genotype plus genotype by environment (GGE) biplots to evaluate genotype performance in specific settings.Genotype and environment vectors are depicted, with a 90 • angle signifying above-average genotype performance in a specific environment.IG5874(G12) underperforms (obtuse angles) across all environments, while IG5863(G8) excels (acute angles) in all but E7.Genotype ranking for E10 places JG14(G42), IG5866(G10), IG5871(G1), IG586(G11), and IG5865(G9) sequentially.Notably, IG5866(G10) and IG5865(G9) maintain yields close to the average, differing from others yielding higher.JG14(G42) shines in E5, whereas ILC239(G40) ranks lowest.Environment performance evaluation (Supplementary Figure S6B) assesses test environments based on IG5871(G1) performance.A line through the biplot origin and genotype marks the axis, depicting environment rankings.IG5871(G1) yields results below average in E10, E3, and E9, near average in E10 and E12, and above average in remaining environments.

Discussion
The integration of analysis of variance and principal component analysis (PCA) has yielded valuable insights into the substantial impacts of location, genotype, and their interactions on yield per plot.Our study underscores the pivotal role of genotype-environment interactions in achieving optimal yields across diverse conditions.The first three PCs explain 86% of the variability, surpassing the minimum reliability threshold of 70% [38,39].In our comprehensive analysis, which accounts for four PCs, the explanatory power extends to 93.2%, demonstrating significant reliability [40][41][42].The AMMI 1 analysis, depicted in the results, provides intricate insights into genotype-environment interactions.This approach facilitates the identification of genotypes excelling in specific environments and those showcasing consistent performance across diverse conditions.Notably, ILC1312(G36) and IG5862(G7) emerge as appealing choices for breeders due to their combined high yield and stability [43][44][45].AMMI 2, a graphical representation using PC1 and PC2 scores, holds distinct advantages over joint-regression-based analysis.It offers a concise overview of complex genotype-environment interactions (GEIs) across multiple settings, elucidating 86% of the G + G × E interaction variation for yield per plot.Thus, the interaction of the 42 genotypes across 12 environments is encapsulated by these three principal components of genotype and environment [46,47].By studying various scenarios, researchers can evaluate how genotype rankings shift with varying weights assigned to productivity and stability.This approach yields valuable insights for genotype selection in breeding programs.To predict yield, a judicious selection between best linear unbiased prediction (BLUP) and additive main effects and multiplicative interaction (AMMI) models is crucial; our findings align with those of Piepho (1994) [48][49][50][51], who concluded that best linear unbiased prediction (BLUP) surpasses any member of the additive main effects and multiplicative interaction (AMMI) family in yield prediction within the context of multi-environment trials (METs) considering intrinsic factors specific to each trial [52][53][54][55].The "which-won-where" pattern suggests diverse cultivar selection according to mega environments.The biplot origin represents an average genotype across environments, with proximity indicating minor G and GE contributions.Longer vectors signify higher G or GE contributions.Noteworthy genotypes with extended vectors (e.g., JG14(G42), IG5866(G10)) exhibit superior performance, while others (e.g., IG5980(G21), ILC1932(G39)) show lower stability [56,57].The angle between genotype vectors and the AEA reflects G and GE components.Certain environments, like E9, E10, and E3, possess notable representativeness and discrimination [58].In this study for genotype selection, we observed associations between test environments, suggesting potential genotype information reduction, leading to cost-effective testing [59].Correlations between closely correlated test environments over time enable the omission of one without significant data loss [60].Our genotype ranking biplot offers critical insights into performance variations across diverse environments, considering mean yield and interactions.This approach holds potential for substantial genotype information reduction, ultimately reducing testing expenses [61,62].

Plant Materials
Forty-one landraces representing the WANA (West Asia and North Africa) region were obtained from the Chickpea Molecular Biology Laboratory, Division of Genetics, IARI, New Delhi (Supplementary Table S7), with one variety, JG14, used as a heat-tolerant check.We cumulatively sowed 160 landraces on the same locations and seasons for heat-tolerance screening (data not presented) and stability analysis, but for better representations of biplots we studied 42 genotypes by using the heat susceptibility index (HSI).The chosen 42 profile genotypes were determined by evaluating grain yield per plot while taking into account the heat susceptibility index (HSI) among all the tested landraces.

Environments and Intercultural Practices
Phenotypic evaluation spanned three locations: ICARDA-Amlaha (23.14711 • N, 76.92035 • E, MSL 502 m) with medium black and alluvial soil, IARI-Regional Research Station-Dharwad (15.45890 • N, 75.00780 • E, MSL 678 m), with tropical black clays derived from the weathering of metavolcanic rocks as block cotton soil, and IARI-New Delhi (28.08 • N, 77.12 • E, MSL 228.61 m), with very deep, somewhat excessively drained alluvial soils.Locations with genotypic code along with latitude, longitude, altitude, year as seasons, date of sowing, and date of harvesting are given in the Supplementary Table S7, with two sowing dates, two seasons (30 days apart), and three locations (2 × 2 × 3) cumulatively representing 12 environments, as coded by E1 to E12.A total of 41 landraces with JG214 check underwent field trials over 2021-2022, including normal and late sowing (30 days apart), and temperature and precipitation of each location and season are given in Supplementary Table S8.An alpha lattice design was used, with three replications per location.Randomized lines/checks occupied 2 m plots.Plants were spaced at 30 cm, rows at 50 cm, plots at 1.5 m, and replications at 2.0 m.Standard agronomic practices were diligently adhered to at each respective location.The background effect of drought was minimized by giving an adequate amount of water by artificial irrigation at all the critical stages of the chickpea crop plant.

Data Collection and Statistical Analysis
The seed yield was recorded after weighing the seeds from a two-meter row plot for each genotype; mean yield with their standard deviations are given in Supplementary Table S9.This study aimed to assess the heat susceptibility index (HSI) and categorize genotypes based on their tolerance to heat stress [37], using heat susceptibility index calculated by (S): (1 − Yh/Yn)/(1 − Xh/Xn), where Yh: yield of genotype under heat stress, Yn: yield of genotype under normal condition, stress intensity = 1 − Xh/Xn, Xh: average yield of all genotypes under heat stress, and Xn: average yield of all genotypes under normal condition.The high-heat-tolerant (HSI < 0.50), heat-tolerant (HSI: 0.51-0.75),moderately heat-tolerant (HSI: 0.76-1.00),and heat-susceptible (HSI > 1.00) scores of all the genotypes were estimated (Fischer and Maurer 1978) using Microsoft Excel.Quantita-tive traits underwent ANOVA to assess variations among genotypes, locations, seasons, and genotype-environment interactions (genotype by location, genotype by season, and genotype by location by season).
The Annicchiarico method, developed in 1992 [63], relies on the recommendation index (ωi(g)), which measures both stability and genotypic adaptability.
Multivariate stability and G × E interaction analysis was performed by using genotype plus genotype by environment (GGE) biplots and AMMI in R Studio (a simplified version of R).Metan facilitated genotype plus genotype by environment (GGE) biplots, while metan managed AMMI analysis [64][65][66], highlighting biplots, and percent of GE interaction was computed from total sum of squares.Genotype plus genotype by environment (GGE) biplots and AMMI graphically showcased G × E interaction and genotype ranking based on mean and stability.Broad-sense heritability was estimated (H 2 = VG/VP) for plot yield using BLUP values across the locations and seasons by dividing genotypic variation to total variation [67].When dealing with high-dimensional data in biplot analysis, both the AMMI and GGE approaches rely on principal component analysis (PCA).However, as the number of components required to capture a significant portion of the original variance increases, it can become challenging to visually interpret the results.In such situations, researchers may need to create multiple biplots to effectively understand the underlying variability in the original genotype-environment interaction (GEI) data (https: //power.larc.nasa.gov/data-access-viewer/(accessed on 10 October 2023)).

Conclusions
Unveiling genotype adaptability and stability using multi-environment trials (METs) of chickpea genotypes provided crucial insights into their adaptability and stability under varying environmental conditions.Notably, the significance of both genotype and interaction effects was remarkably pronounced (p < 0.001).In particular, the evaluation of genotypes ILC239(G40), IG6002(G26), IG5862(G7), ILC0(Czech)(G30), and IG5861(G6) unveiled their remarkable adaptability, with the lowest weighted average absolute scores of BLUP (WAASB) values (0.77, 0.93, 1.03, 1.10, 1.21) observed when considering four or more interaction principal component axes (IPCAs).These genotypes demonstrated resilience across a spectrum of scenarios, as affirmed by the multifaceted statistical results derived from genotype plus genotype by environment (GGE) and AMMI biplots.Comparatively, best linear unbiased prediction (BLUP) outperformed any member of the AMMI family in accuracy.The multi-environment trials (METs) approach unraveled the intricate dynamics of chickpea genotypes, illuminating their versatility and robustness.This study lays a solid foundation for informed genotype selection and breeding advancements, aiming towards bolstered chickpea yields and resilience in the face of varying environmental challenges.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/plants12213691/s1, Figure S1 S1: List of selected genotypes used in this study with their heat susceptibility index values of individual environment; Table S2: Genotype rankings in respective tested environments; Table S3: Annichiarico environmental index; Table S4: Descriptive statistics, variance components, and genetic parameters for grain yield of 42 genotypes evaluated across the year and locations; Table S5: Genotype code and their yield across environments with their IPCA scores; Table S6: Genotypes mean yield with different stability indices for ranking genotypes under multi-environmental condition; Table S7: Environmental description of the experimental sites; Table S8: Mean monthly temperature and precipitation data of all the locations; Table S9: Mean ± standard deviation of yield traits of individual locations and seasons.

Figure 1 .
Figure 1.Violin plots of yield trait of each individual environment are distributed for both normal (blue color) and late (green color) conditions.

Figure 1 .
Figure 1.Violin plots of yield trait of each individual environment are distributed for both normal (blue color) and late (green color) conditions.

Figure 3 .
Figure 3. (A) The biplot presents the relationship between grain yield and the weighted average of absolute scores for the best linear unbiased prediction of genotype-environment interaction (WAASB) for 42 genotypes tested across 12 environments.A black circle in the IV quadrant represents a theoretical genotype with both high productivity and broad adaptability.Horizontal and vertical black arrows indicate the directions of yield increase and stability improvement, respectively.(B,C) The comparative predictive precision of the AMMI (additive main effects and multiplicative interaction) family and best linear unbiased prediction (BLUP) methods for the yield of 42 C. arietinum genotypes are depicted.The distribution of 1000 estimates of root mean square prediction difference (RMSPD) is visualized through boxplots.

Figure 3 .
Figure 3. (A) The biplot presents the relationship between grain yield and the weighted average of absolute scores for the best linear unbiased prediction of genotype-environment interaction (WAASB) for 42 genotypes tested across 12 environments.A black circle in the IV quadrant represents a theoretical genotype with both high productivity and broad adaptability.Horizontal and vertical black arrows indicate the directions of yield increase and stability improvement, respectively.(B,C) The comparative predictive precision of the AMMI (additive main effects and multiplicative interaction) family and best linear unbiased prediction (BLUP) methods for the yield of 42 C. arietinum genotypes are depicted.The distribution of 1000 estimates of root mean square prediction difference (RMSPD) is visualized through boxplots.

Plants 2023 , 17 Figure 4 .
Figure 4. (A) The heatmap illustrates the rankings of 42 chickpea genotypes based on the utilization of interaction principal component axes (IPCAs) in the weighted average of absolute scores for the best linear unbiased predictions (BLUPs) of genotype-environment interaction (WAASB) estimation.(B) Circos plot showing genotypes contribution to individual IPCAs.(C) Circos plot of individual genotypes contribution to yield and respective IPCAs.

Figure 4 .
Figure 4. (A) The heatmap illustrates the rankings of 42 chickpea genotypes based on the utilization of interaction principal component axes (IPCAs) in the weighted average of absolute scores for the best linear unbiased predictions (BLUPs) of genotype-environment interaction (WAASB) estimation.(B) Circos plot showing genotypes contribution to individual IPCAs.(C) Circos plot of individual genotypes contribution to yield and respective IPCAs.

Plants 2023 , 17 Figure 5 .
Figure 5. (A) The environment vector perspective of the GGE biplot illustrates the resemblances among test environments.(B)The GGE biplot polygon view depicting the "which-won-where" pattern showcases the combined impact of genotype main effects and G × E interaction for 42 C. arietinum genotypes across two seasons and three locations in terms of yield per plot.These biplots were generated using centering = 2, SVP = 3, and scaling = 0 settings.

Figure 5 .
Figure 5. (A) The environment vector perspective of the GGE biplot illustrates the resemblances among test environments.(B)The GGE biplot polygon view depicting the "which-won-where" pattern showcases the combined impact of genotype main effects and G × E interaction for 42 C. arietinum genotypes across two seasons and three locations in terms of yield per plot.These biplots were generated using centering = 2, SVP = 3, and scaling = 0 settings.

Figure 6 .
Figure 6.(A) The GGE biplot s depiction of the "mean vs. stability" pattern.(B) The GGE biplot "genotypes ranking" pattern for genotype comparison with ideal genotype showing G + G × E interaction effect of 42 genotypes under two seasons and three locations for yield per plot.The biplots were created based on centering = 2, SVP = 2, and scaling = 0.

Figure 6 .
Figure 6.(A) The GGE biplot's depiction of the "mean vs. stability" pattern.(B) The GGE biplot "genotypes ranking" pattern for genotype comparison with ideal genotype showing G + G × E interaction effect of 42 genotypes under two seasons and three locations for yield per plot.The biplots were created based on centering = 2, SVP = 2, and scaling = 0.

Figure 7 .
Figure 7. (A) GGE biplot "discriminativeness vs. representativeness" pattern compares genotypes, highlighting the ideal genotype with G + G × E interaction.(B) GGE biplot "Env.Ranking" pattern compares environments, highlighting ideal environment with G + G × E interaction for 42 genotypes, 2 seasons, and 3 locations, using yield per plot.Biplots: centering = 2, SVP = 2, and scaling = 0. 2.6.5.Test Environment Discrimination and Representation Figure 7B illustrates distances among test environments for all traits.Concentric circles in the biplot visualize environment vector lengths, proportional to standard deviations within environments, reflecting environment discrimination capability.E5 and E6 are most informative, whereas E10 and E9 are least discriminative.Representativeness of test environments: E10 is highly representative, unlike E5 and E6 which are less so.Environments that are both discriminative and representative, like E10, suit broad genotype selection.Nonrepresentative yet discriminative test environments (e.g., E5, E6) aid specific megaenvironment genotype selection.From Figure 7B rankings, preferred target

:
Plots showing the pictorial representation of genotype and overall mean performance of studied traits at each environment under study; Figure S2: (A) AMMI 2 biplots of PC1, PC2, and PC3; Figure S3: Calculated values of weighted average of stability (WAASB) and mean performance (Y) (WAASBY); Figure S4: (A) Heatmap of genotypes × environmental interactions, (B) correlation between stability indices; Figure S5: (A) Comparative analysis of JG14 (G42) and ILC1932 (G39) across diverse environments (SVP = 3), (B) genotype plus genotype × environment (GGE) biplot "environment × genotype relationship" view; Figure S6: (A) Genotype ranking according to performance in a specific environment (E10), (B) ranking trial environments based on the relative performance of genotypes G1 and G2; Figure S7: Field pictures of the heat tolerant testing; Figure S8: Graph of monthly average temperature of chickpea growing season of 2021 to 2023.Table

Table 1 .
Estimation of significant level for yield and yield contributed traits of 42 C. arietinum accessions revealed by ANOVA.

Table 1 .
Estimation of significant level for yield and yield contributed traits of 42 C. arietinum accessions revealed by ANOVA.